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Globally coupled phase oscillator models, such as the Kuramoto model, exhibit spontaneous col¬ 
lective synchronization. Such models can be restated in terms of interactions within and between 
subsets of oscillators. An approximation for the internal structure of coherent subsets of oscillators 
can be made based on the observations of natural ordering and uniform tidal effects. The approx¬ 
imation is seen to perform well for predicting the microstructure in a variety of phase oscillator 
models. 


Spontaneous collective synchronization among nearly 
identical phase oscillators is a classic emergent phe¬ 
nomenon and phase transition [l|. The first theoretical 
model to provide insight into this transition was provided 
by Winfree in 1967 [3 and takes the form 

K ^ 

1=1 

The functions P and Q must be 27r-periodic, and so can 
be expressed as trigonometric functions of angular differ¬ 
ences and sums. Kuramoto considered a simpler coupling 
m given by: 

K ^ 
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The simplest coupling is T (A0) = sin (A0). Such a cou¬ 
pling tends to draw two nearby oscillators together, while 
the disorder in the natural rotation rates, LOi, tends to 
spread them out. For many choices of P and Q or of T, 
if the coupling K is strong enough and the variation of 
the uji is small enough, a subset of the population will mu¬ 
tually entrain: they will lock into relative positions that 
are essentially hxed after an initial transient. Winfree 
argued that the fraction of mutually entrained oscillators 
exhibits a phase transition in the coupling K, and Ku¬ 
ramoto analytically calculated the magnitude of the en¬ 
trained fraction for his model. Various models have been 
considered by other authors and together they have 
been used to describe the collective behavior of many 
physical systems including coupled Josephson junctions 
Q, laser arrays Q, electronic auto-oscillators Q, elec¬ 
trochemical oscillators [1, 0, and acoustically coupled 
mechanical rotors (Til . [T3 . 

Globally coupled phase oscillator models can be rewrit¬ 
ten exactly as mean field models. Because of this, most 
research has focused on approximating the behavior of 
the mean field. The most important results are found 
with the approximation of large system size, V —>■ 00 
fllillil. Work focusing on finite populations relies 
on ensemble averaging to make any meaningfully general 
claims [iM3- Unfortunately, the only published mech¬ 
anism for elucidating population-specific tipping points 
is to measure or simulate the behavior directly. 


Figure 1: Pairwise correlations between all oscillators in a 
simulation of the Kuramoto model, N = 64. Lighter and 
darker intensity indicates smaller and larger values of pij, re¬ 
spectively. The indices were not optimized for block diagonal 
structure but were simply sorted according to uji. For this 
simulation, the coupling strength is equal to the width of the 
population’s Gaussian distribution, placing it well below the 
critical coupling. 

Presence of coherent subsets The mean field reformu¬ 
lation is exact, but it theoretically obscures an important 
observation. Figure [1] illustrates how groups of oscillators 
clump into groups (such as explored by Maistrenko et al. 
0)- The intensity indicates the value of pij , defined as 



using i = V—1- The obvious block diagonal structure is 
obtained simply by sorting the indices according to natu¬ 
ral frequency uii. Subsets of oscillators with nearby natu¬ 
ral frequencies mutually entrain. Although not indicated 
in this figure, oscillators in these clumps also appear to 
act in concert, suggesting that their degrees of freedom 
are not independent. Taken together, these observations 
indicate that finite-sized phase oscillator models exhibit 
rich internal structure. 
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Figure 2: Depiction of two subsets (black dots) in a full pop¬ 
ulation (black and gray dots) of 64 oscillators, not related to 
the simulation in figure [T] Subset 7 is more spread out than 
subset 4, and so has a smaller mean field relative to subset 
size. If these subsets are coherent, their constituent oscilla¬ 
tors may get further apart or closer together, but they do not 
lap each other. 


The presence of coherent subsets gives rise to a num¬ 
ber of interesting questions. How do we best quantify 
the behavior of these coherent subsets? What approx¬ 
imations must be made to express the dynamics of the 
model in terms of the subsets’ dynamics? How do the 
subsets evolve? Can we approximate the structure or 
evolution of the subsets? All of these questions will be 
addressed in this letter. 

Subset reformulation The first papers on synchro¬ 
nization established that globally coupled phase oscilla¬ 
tor models can be written exactly in terms of mean fields 
ii. For models with coupling that depends only on 
the first harmonic, such as the Kuramoto model, we use 
the mean field 
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(Models that use the harmonic need to consider the 
mean field defined using By a trick of complex 

algebraflj, the dynamics of an individual oscillator can 
be exactly rewritten in terms of interactions with the 
mean field. For the Kuramoto model, the interactions 
take the form 


9i = uji +r K - 6i). (5) 

If one could predict the behavior of r and ip, one would 
essentially solve the dynamics of the Kuramoto model 
and could predict the properties of interest for the sys¬ 
tem. Ott and Antonsen provide a method for perform¬ 
ing such a calculation for the Kuramoto model under the 
assumption of a continuous oscillator density M- Un¬ 
fortunately, there is no generic closed form solution for 
r (t) or Ip (t) for an arbitrary finite population behaving 
according to the Kuramoto model or other similar mod¬ 
els. 

How might we improve upon the mean-field reformu¬ 
lation for finite sized systems? Consider the following 


minor refinement: split the summation of the traditional 
order parameter into Af pieces: 


N 
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In this expression, the oscillator’s contribution 
occurs only once on the left and right; that is, each oscil¬ 
lator is assigned to only one subset, S£. The assignment 
of oscillator di to subset si is (for the moment) arbitrary 
PH . Taking a cue from the traditional analysis, define 
the subset mean field for the subset as 

^ (7) 

jest 



Two such subsets, along with their mean fields, are de¬ 
picted in figure Employing the geometric interpreta¬ 
tion of complex numbers, we see that ipi points roughly 
to their average phase. If a subset’s oscillators are near 
to each other then Rg Ri Ng and if they are scattered 
about the unit circle then Rg <C Ng. The subset’s mean 
field serves as a simple measure of the coherence and lo¬ 
cation of the subset, and its formulation emerges from 
the original order parameter definition. 

The traditional order parameter can be calculated, 
without any approximations, in terms of these mean 
fields. A simple algebraic substitution leads to 
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The behavior of each oscillator can also be written ex¬ 
actly in terms of interactions with these mean fields. The 
specifics depend upon the functions P and Q in equation 
[T] or r in equation but such a restatement can be ob¬ 
tained in terms of the harmonic mean fields after employ¬ 
ing the proper trigonometric identities. For example, the 
oscillator dynamics for the Kuramoto model are 

K ^ 

0i=^i +Re sin {ipg-0,), (9) 

^ e=i 


which closely resembles the mean-field interaction given 
in equation [SI Again, this equation does not rely on any 
approximations. The only limitation is that I have yet to 
predict Rg and ipg. 

In order to predict the behavior of the subset mean 
fields, we need an expression for their dynamics. Eval¬ 
uating the time derivative of the definition (equation [7]) 
and rearranging leads to 


Rg + iTpgRg = (10) 

jest 

= Oj sin {ipg — 9j) + 1 Y^ ^j cos {ipg — Oj). 

jest jest 

( 11 ) 
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The real and imaginary components give separate equa¬ 
tions for the evolution of a subset amplitude and phase, 
respectively. To obtain the equations for a specific model, 
Oj must be replaced with the expression for the time 
derivative specific to that model. Although the dynam¬ 
ics in 9j depend upon all of the mean fields, judicious use 
of trigonometric identities can lead to separable contribu¬ 
tions for internal dynamics and interactions between sub¬ 
set mean fields. The internal dynamics are of particular 
interest: if these could be approximated in terms of the 
subset’s mean field, it would be possible to coarse-grain 
the dynamics of the oscillator model into a model for in¬ 
teracting amplitude oscillators representing the subsets. 
Such a coarse graining relies upon a good approximation 
for tpi — 9j. 

Coherent subset approximation One good approxima¬ 
tion for ipi — 9j that is independent of model details is 
the coherent subset approximation. The approximation 
is built on two observations about the clumps discussed 
on the first page: 

1. Natural ordering: For oscillators in a coherent sub¬ 
set, the positions 9i order according to increasing 
natural speed, uji. 

2. Uniform tidal effects: Relative oscillator positions 
fluctuate nearly in unison. For oscillators in a co¬ 
herent subset £, relative positions can be related to 
the subset’s mean field amplitude. Re. 

While uniform tidal effects have not previously been 
noted, evidence of natural ordering is abundant in the 
literature, beginning with Winfree’s discussion of synta- 
lansis in his original paper on the topic j^. One math¬ 
ematical approximation that reflects these two observa¬ 
tions is 

9, (ujj - wef (W - Rtf ■ (12) 

The w-dependence reflects the observation of natural or¬ 
dering. The dependence on the subset’s mean field am¬ 
plitude and the neglect of other time dependence reflects 
the observation of uniform tidal effects. 

In the simplest calculation, the constant C and expo¬ 
nents OL and /3 do not depend on the details of the phase 
oscillator model, but only on the definition of the sub¬ 
set mean field. With a small rearrangement, equation [3 
becomes 

^ (13) 

To obtain constraints on the constant and exponents, ex¬ 
pand the real and imaginary components of equation 1131 
to second order in 9j — ipt- The imaginary component 
gives 

0 = ^ sin {9j - ipt) 

^C{Nt- Rtf Y. ’ (14) 



Predicted 9 (rad) 

Figure 3: Histogram of actual (y) vs predicted (x) relative 
phase positions for the approximation given by equations [16] 
for the Kuramoto model. Intensity is linear in density. Figure 
(a) shows results for small subsets, each of which constitute 
fewer than 10% of their simulation’s population. Figure (b) 
shows results for large subsets, each of which have more than 
10% of their simulation’s population. The data represent ag¬ 
gregated results of many simulations with a variety of coupling 
strengths and full population sizes, and with emergent sub¬ 
sets of various sizes. The correlation between the prediction 
and the measurement is — 0.86 for (a), = 0.99 for (b). 

which can only hold generically if a = 1. The real com¬ 
ponent gives 

Rt = Y (^J “ 

jese 

^Ar,-^(7V,-R,)2^^(a;,-a;,)"“, (15) 

jesi 

which implies ,5 = 1/2 and C~‘^ = ^ J2jesi (^1 ~ — 

Aj. Finally we arrive at 

9j -i>t- ^Nt - Rt, (16) 

At 

This equation is the key result of this letter. Notice that 
the derivation does not rely on a particular form of P, 
Q, or F. Although this equation does not predict which 
oscillators will aggregate into coherent subsets, it gives 
a mo del-independent time-dependent prediction for the 
steady-state structure of such subsets. 

Agreement with simulation The agreement between 
the actual relative phases for the Kuramoto model and 
the predictions of equation |TH| is shown in the density 
plots given in figure [3] and its inset. The main figure 
represents data from small subsets while the inset repre¬ 
sents data from large subsets. The figure includes data 
from many different empirically identified subsets, them- 









4 


selves constituents of about 80 unrelated unimodal pop¬ 
ulations of various sizes, and all simulated with five cou¬ 
pling strengths. The coupling strengths are chosen ran¬ 
domly between 0 and 2; the critical coupling strength for 
these simulations is roughly 1.6 [l|. Additional details as 
well as other phase oscillator models are discussed in the 
supplementary material. 

Figure [3] illustrates a strong correlation between the 
measurement and prediction: the approximation works 
well. Other models, such as the Ariaratnam-Strogatz 
model 01 or a third harmonic approximation to the saw¬ 
tooth function, exhibit equally impressive small-subset 
results. The prediction is not perfect, with nonlineari¬ 
ties particularly prominent in the inset at large relative 
positions. It is perhaps unsurprising that large-subset 
data exhibit model-specific nonlinearities: these agree 
with the predictions of Sakaguchi and Kuramoto 01 . It is 
remarkable, however, that the nonlinearities are one-to- 
one. This one-to-one nature strongly suggests that more 
accurate predictions can be expressed as model-specific 
power series of equation 1161 In other words, for models 
that exhibit natural ordering and uniform tidal effects, 
equation 1161 serves as a universal starting point. 

Discussion and conclusion The coherent subset ap¬ 
proximation differs from similar results in the literature 
in key ways. Sakaguchi and Kuramoto predicted en¬ 
trained relative phases for their model in ref 0| equa¬ 
tion 5a. Their predictions only applied to oscillators 
entrained to the order parameter whereas equation [13] 
describes all coherent behavior, whether the mean field 
is appreciable or negligible, and whether or not the os¬ 
cillators of interest are entrained to the mean field. It 
also gives a prediction for the instantaneous relative po¬ 
sition which does not rely on long-time averaging or the 
coupling strength. Others have predicted the oscillator 
density distribution, p(0,a;,t), for the Kuramoto model 
with infinite system size and various population distri¬ 
butions [1^ 13, . Equation [131 is distinct from those 

predictions by addressing finite-size synchronization for 
a broad set of models, and by imposing no restriction 
on the population distribution. The coherent subset ap¬ 
proximation is both more precise and more general than 
previous predictions. 

The coherent subset approximation provides a new an¬ 
gle for analyzing distributions with finite support. In 
their impressive analysis, Martens et al. analytically 
compute nearly all of the bifurcation diagram for the Ku¬ 


ramoto model with a Lorentzian bimodal frequency dis¬ 
tribution (^ . As explained in the supplementary mate¬ 
rial, obtaining the saddle-node infinite-period (SNIPER) 
bifurcation curve for small a is simple using the coherent 
subset approximation, and leads to unexpected scaling. 
Martens predicted that ct oc 2 — wq as wq —t 2. The co¬ 
herent subset approximation, on the other hand, predicts 
that a oc y/2 — dig as wq 2. The range of values over 
which the scaling takes this form is not yet clear, but it 
should hold for any sufficiently narrow distribution with 
finite support. Whether an infinite bimodal Gaussian 
distribution would have linear or square-root scaling is 
intriguing but unknown. 

All globally coupled phase oscillator models can be de¬ 
composed. The evolution of the subset mean fields is 
given by equation [TT] even those models which do not 
exhibit natural ordering and uniform tidal effects. Un¬ 
fortunately, the mean field dynamics given by equation 
[TT] depend on model-specific details and so do not coarse 
grain to a universal form. If universality exists across a 
broad set of phase oscillator models, it will arise because 
the statistics of the aggregations do not depend upon the 
underlying model. This work provides the basis for such 
an analysis, but any such claims go beyond the scope of 
this letter. 

In this letter, I have shown that the dynamics of glob¬ 
ally coupled phase oscillator models can be restated ex¬ 
actly in terms of interactions among and within subsets. 
Phase oscillator models that exhibit natural ordering and 
uniform tidal effects will exhibit coherent subsets they 
have nearly identical microstructure. The approximation 
given in equation [T3] appears to perform well and could 
serve as a starting point for more accurate approxima¬ 
tions in a wide variety of models. Rather than relying 
on long-time or ensemble averaging to make meaning¬ 
ful predictions, an analysis based on subsets provides a 
formulation for predicting features of individual popu¬ 
lations, paving the way for insights into a decades old 
problem in dynamic phase transitions. 
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I. NUMERICAL SETUP 


For all simulations, the population is sampled from a normal distribution of width a = 
and mean u = The mean is inconsequential for the Kuramoto model but 

matters for one of the models considered below. In all simulations, I choose positions that 
are randomly distributed around the unit circle. The number of oscillators for a given 
simulation is sampled randomly from the uniform distribution between 30 and 5000. 

I choose K values randomly from the uniform distribution between 0 and 2. Other models 
use other ranges as detailed below. The coupling ranges are based on trial and error, aiming 
to get a good sampling of unsynchronized and synchronized behavior for each model. 

In all simulations presented, I follow the well-established practice of using a fourth-order 
Runge-Kutta method with a time step of 0.01s til. Based on the numerical evidence of 

A 

Chulho et ah (l|, the transient behavior of the Kuramoto model lasts for about 
seconds. I run all simulations for 900\/]V^ seconds, i.e. 30 multiples of the transient 
duration, storing the order parameter and oscillator positions each second, i.e. every 100 
time steps. I then use the Mean Squared Error Rule to identify when the order parameter 
reaches steady state Q, using either that or 3y/Nose seconds, whichever is larger. 

II. IDENTIFYING ENTRAINED SUBSETS OF OSCILLATORS 

In addition to the order parameter and oscillator positions, I also monitor the phase 
differences between oscillators with consecutive natural speeds, and make note of all 27r 
phase slips and the time step at which they occur. To identify coherent subsets, discard all 
phase slips that occur before the onset of the steady state and identify all pairs with no slips 
for the duration of the steady state. Such pairs are assumed to be mutually entrained. If 
one oscillator is entrained to two other oscillators, then all three of them must be mutually 
entrained. Using this approach, I extract the entrained subsets empirically. 

III. PRODUCING THE THIRD FIGURE 

To produce the third figure in the letter, I select 1000 random moments in time, 
during the steady-state of each subset. I compute the subset mean field, (f„) 
from the 9j (t„). I then compute the actual relative position, 9j {tn) — ‘ipe {tn), as well as 
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Figure 1: Prediction vs simulation for the Kuramoto model. Figure (a) is for small coherent subsets 
while figure (b) is for large coherent subsets. In contrast to the figure in the letter, the intensity 
indicates the logarithm of the density, rather than the density itself. 



Predicted 6 (rad) Predicted 9 (rad) 


Figure 2: Prediction vs simulation for the Ariaratnam-Strogatz model, a Winfree model. Figure 
(a) is for small coherent subsets while figure (b) is for large coherent subsets. The correlation for 
the small subsets is = 0.99 and for the large subsets is = 0.99. 


the prediction for the relative positions based on The plots are two-dimensional 

histograms depicting the number of times a predicted value of 9j — ijj£ for any one oscillator 
corresponded with any measured value of 6j — ipe. 


IV. ADDITIONAL MODELS 

In this section I present the results of the Kuramoto model alongside the results of two 
other globally coupled phase oscillator models. 
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Predicted 6 (rad) Predicted 6 (rad) 

Figure 3: Prediction vs simulation for the third-harmonic sawtooth model, a Kuramoto-like model. 
Figure (a) is for small coherent subsets while figure (b) is for large coherent subsets. The correlation 
for the main figure is = 0.86 and for the inset is = 0.95. 

Figured] shows the small and large histograms for the Kuramoto model, using the same 
data as in the letter. In contrast to the hgures in the main letter, the intensity indicates 
the logarithm of the density. This is necessary to represent the results of the sawtooth 
model, discussed below. Even with heightened sensitivity to small densities, the interpre¬ 
tation remains the same: relative positions are well approximated by the coherent subset 
approximation for small subsets, and one-to-one nonlinearities arise for large subsets. 

Ariaratnam and Strogatz considered a Winfree model for which P {6j) = —1 — cosOj 
and Q {6i) = sin6'i. Figure [2] shows the small and large subset histograms for the model. For 
these simulations, the coupling was chosen randomly from the uniform distribution between 
0 and 1. Strong correlations and nonlinear effects occur for subsets greater than 3% of the 
full population size, substantially smaller than for the Kuramoto model. 

Note that the mean frequency is signihcant in this model. This model can be interpreted 
as a Kuramoto-like model in which the oscillators feel a pull between a “uniform held” 
directed toward 6 = 0 coupled with strength K, a mean held coupled with strength iF/2, 
and an odd additive term—sin (0* + —^Iso coupled with strength iF/2. We expect the 

mean held to have two elements: a component aligned to the uniform held, and a component 
rotating in spite of the uniform held. Large u should lead to more oscillators in the rotating 
component. Ariaratnam and Strogatz considered an evenly spaced distribution of oscillators, 
without disorder, so it is not known exactly how the choice of u will ehect the results, but 
an ehect is expected. 
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The final model is a Kuramoto-like model based on the Fourier expansion of the sawtooth 
function truncated at the third harmonic, for which F (A0) = sin A6 — ^ sin 2A6 + | sin 3A6. 
The results for this model are shown in figure [3l The coupling was chosen randomly from 
the uniform distribution between 0 and 4. The cutoff for large subsets was empirically set 
at 45%, substantially larger than for the Kuramoto model. It is also notable that the vast 
majority of the small-subset density lies within very small relative angles. 

As mentioned in the main text, the predictions work fairly well, and although nonlineari¬ 
ties are present, the predictions still seem to have a one-to-one mapping to the measurements 
for most data. It is particularly interesting to note that large subset nonlinearities for the 
Kuramoto and the third-order sawtooth closely resemble the actual form of the coupling 
functions themselves! 


V. SQUARE-ROOT SCALING OF THE SNIPER BIFURCATION CURVE 


In this section I will obtain a finite-size correction to the phase diagram of Martens et 


al 


Q. 


To begin, I insert the form of 6j from the Kuramoto model into equation 11 in the 
letter. I also substitute the coherent subset approximation, equation 16 in the letter, and 
expand the trigonometric sums to second order to obtain 


Re = -2 Ae\/Ne - Re + 2{Ne- Re) cos (-0™ - '<Pi ), (1) 

m 

K X 

■ipeRi = AeRe + {2 Re- Ne) sin {ijm - i’e) ■ (2) 

m 

For a symmetric bimodal population, for which all the oscillators are members of one of 
two coherent subsets Si or S 2 ; I can obtain an expression for the mean field amplitudes 
Ri = R 2 = Re as well as the difference of the phases, A'lp = '<p 2 — '>p 2 '- 


Re — —2Ae\/Ne — Re + 2 (iV^ — Re) j^R^ A cos At):) , 

^ A / - - r,‘^ R^ ~ L- ■ A / 

—A'tp = ijj 2 — oJi — 2 -—-A sm A'lp. 

at N 


(3) 

(4) 


Now suppose that the peaks of the bimodal distribution are very narrow with respect to 
the coupling strength and their separation. In that case, the oscillators belonging to each 
peak should mutually entrain, behaving like two giant oscillators. The saddle-node infinite- 
period (SNIPER) bifurcation separates the the phase space where these two giant oscillators 
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are locked together or are counter-propagating. The dynamics for the locked oscillators are 
simply the hxed points given by = 0 and = 0, so 

2A£a /— Ri = 2 (Nf — Ri) (1 + cos A'^), (5) 

0)2 — a)i = 2- ^ —^ifsinA'^. (6) 

The saddle-node boundary occurs at values of A^ and 002 — oJi at which these equations can 
not be satished. The most extreme separation between 0)2 — ui occurs when sin A'^ 1, in 

which case cos A'^ 0, leading to 

2A£a /— Ri = 2 {Ni — Re) (7) 

U2-uJi = 2- — - K. (8) 

To obtain expressions using the notation of Martens, use 


0^2 — (Ui — 2 CJq, 


a2 _ 1 2 

Af — (^e- 


Substituting this into the above expressions and eliminating Re leads to 

Aae 


__ 1 / 9 _ 

K K 


2 + 


Auq 


2N 

N-2' 


Taking the large N limit so that N/ [N — 2) ^ 1, and using Martens’ dehnitions, 


a = 

Aae 

K '' 


4 Wo 

cDo = 

K 


(9) 

( 10 ) 


( 11 ) 


( 12 ) 

(13) 


I obtain 

/x 

d = a/ 2 — (ho (2 -|- (ho). (14) 

In contrast, in the vicinity of Cjq ~ 2, Martens’ eq 33 comes to 

A = ^(2-(ho), (15) 

where A represents the half-width-half-maximum of the individual peaks and A = 4 A/i^. 

I observed in the previous section that the coherent subset approximation is not exact 
for large separations 'ipe ~ dj. These inconsistencies limit the range of applicability of my 
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calculation, but do not invalidate it. For subsets with sufficiently narrow d, all relative 
positions should fall within the well-ht range of the prediction. The square-root dependence 
will certainly hold for these distributions. 
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